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Abstract 

Multitype branching processes (MTBP) model branching structures, where the nodes 
of the resulting tree are objects of different types. One field of application of such models 
in biology is in studies of cell proliferation. A sampling scheme that appears frequently is 
observing the cell count in several independent colonies at discrete time points (sometimes 
only one). Thus, the process is not observable in the sense of the whole tree, but only 
as the "generation" at given moment in time, which consist of the number of cells of 
every type. This requires an EM-type algorithm to obtain a maximum likelihood (ML) 
estimation of the parameters of the branching process. A computational approach for 
obtaining such estimation of the offspring distribution is presented in the class of Markov 
branching processes with terminal types. 
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1 Introduction 

Branching processes (BP) are stochastic models used in population dynamics. The theory of 
such processes could be found in a number of books ( [12] , [2] , pQ ) , and application of branching 
processes in biology is discussed in [H], [22], [15], [TT]. Statistical inference in BP depends 
on the kind of observation available, whether the whole family tree has been observed, or 
only the generation sizes at given moments. Some estimators considering different sampling 
schemes could be found in [10] and [23]. The problems get more complicated for multitype 
branching procesess (MTPB) where the particles are of different types (see [IE])- Yakovlev 
and Yanev |23] develop some statistical methods to obtain ML estimators for the offspring 
characteristics, based on observation on the relative frequencies of types at time t. Other 
approaches use simulation and Monte Carlo methods (see [E], [8]). 

When the entire tree was not observed, but only the objects existing at given moment, an 
Expectation Maximization (EM) algorithm could be used, regarding the tree as the hidden 
data. Guttorp |10| presents an EM algorithm for the single-type process knowing generation 
sizes and in [9] an EM algorithm is used for parametric estimation in a model of Y-linked gene 
in bisexual BP. Such algorithms exist for strictures, called Stochastic Context-free Grammars 



(SCFG). A number of sources point out the relation between MTBPs and SCFGs (see [2"U] . 
[7]). This relation has been used in previous work [5] to propose a computational scheme for 
estimating the offspring distribution of MTBP using the Inside-Outside algorithm for SCFG 
(|16]). A new method, related to this, but constructed entirely for BP will be presented here. 

The EM algorithm specifies a sequence that is guaranteed to converge to the ML estimator 
under certain regularity conditions. The idea is to replace one difficult (sometimes impossible) 
maximization of the likelihood with a sequence of simpler maximization problems whose limit 
is the desired estimator. To define an EM algorithm two different likelihoods are considered - 
for the "incomplete-data problem" and for "complete-data problem". When the incomplete- 
data likelihood is difficult to work with, the complete-data could be used in order to solve 
the problem. The conditions for convergence of the EM sequence to the incomplete-data ML 
estimator are known and should be considered when such an algorithm is designed. More 
about the theory and applications of the EM algorithm could be found in [17j . 

The paper is organized as follows. In Section 2 the model of MTBP with terminal types 
is introduced. Section 3 shows the derivation of an EM algorithm for estimating the offspring 
probabilities in general, and then proposes a recurrence scheme that could be used to ease the 
computations. A comprehensive example is given in Section 4 and the results of a simulation 
study are shown in Section 5. 

2 The Model 

A MTBP could be represented as Z(i) = (Z x (t), Z 2 (t), . . . Z d (t)), where Z k (t) denotes the 
number of objects of type T& at time t, k = 1, 2, . . . d. An individual of type k has offspring of 
different types according to a d-variate distribution pk(xi, x 2 , ■ ■ ■ , Xd) and every object evolves 
independently. If t = 0,1,2,..., this is the Bienayme-Galton- Watson (BGW) process. For 
the process with continuous time t G [0, oo), define the embedded generation process as follows 
(see [2]). Let Y n = number of objects in the n-th generation of Z(t). If we take the sample 
tree 7r and transform it to a tree n' having all its branches of unit length but otherwise 
identical to n, then Y n (7r) = Z n (7r'), where Z n is a BGW process. We call Y n the embedded 
generation process for Z(i). The trees associated either with BGW process, or the embedded 
generation BGW process will be used to estimate the offspring probabilities. 

Now we consider MTBP where certain terminal types of objects, once created, neither 
die nor reproduce (see [20]). If T = {T\,T 2 , . . . ,T m } is the set of non-terminal types and 
T T = {T^ ,Tj ', . . . ,Tj_ m } is the set of terminal types, then an object of type Tj produces 
offspring of any type and an object of type Tj does not reproduce any more. Here each 
Tj G T t constitutes a final group (see p2]). This way we model a situation where "dead" 
objects do not disappear, but are registered and present as "dead" through the succeeding 
generations. The process described above is reducible, because once transformed into a 
terminal type, an object remains in it's final group. For irreducible processes some statistical 
theory has been developed (see [3] for example), but for reducible ones statistical inference 
is case-dependent. 

We are interested in estimation of the offspring probabilities. If the whole tree ir is 
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Figure 1: A discrete time process a), continuous time process b) and its embedded process 
c). 



observed the natural ML estimator for the offspring probabilities is 



p(T v -+A) = 



c(T v -» A) 
c(T v ) ■ 



(1) 



where c(T v ) is the number of times a node of type T v appears in the tree ir and c(T v — > A) 
is the number of times a node of type T v produces offspring A in ir. It is not always possible 
to observe the whole tree though, often we have the following sampling scheme {Z(0), Z(t)}, 
for some t > 0. Let Z(0) consists of 1 object of some type. Suppose we are able to observe a 
number of independent instances of the process, meaning that they start with identical objects 
and reproduce according to the same offspring distribution. Such observational scheme is 
common in biological laboratory experiments. If t is discrete Z(i) is the number of objects 
in the t-th. generation. For continuous time it is a "generation" in the embedded generation 
BGW process. A simple illustration is presented in fig. 1 a)— c), where "alive" objects are 
grey, "dead" ones are white and numbers represent the different types. 

Here the notion that "dead" objects present and could be observed in succeeding gen- 
erations as terminal types is crucial. If "dead" particles disappeared somewhere inside the 
"hidden" tree structure, estimation would be impossible (see fig. 2 for an example). 
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Figure 2: Because of disappearing of a particle of type 2 in the first tree information about 
the reproduction has been lost and we have the same observation as in the second tree. 
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3 The EM Proposal 



3.1 The EM Algorithm 

The EM algorithm was explained and given its name in a paper by Dempster, Laird, and 
Rubin in 1977 [6J. It is a method for finding maximum likelihood estimates of parameters in 
statistical models, where the model depends on unobserved latent variables. Let a statistical 
model be determined by parameters 9, x is the observation and Y is some "hidden" data, 
which determines the probability distribution of x. Then the density of the "complete" 
observation is f(x,y\9) and the density of the "incomplete" observation is the marginal one 
g(x\8) = J f(x,y\9)dy. We can write the likelihoods and the conditional density of Y given 
x and 9 this way: 

f(r t\9) 

L(9\x,y) = f(x,y\8), L(9\x) = g(x\9), and k(y\9,x) = JK ' . 

g(x\9) 

The aim is to maximize the log likelihood 

\ogL{9\x) = logg(x\9) = log L(0\x,y) - log k(y\9,x). 

As y is not observed the right side is replaced with its expectation under k(y\9',x): 

log L{8\x) = Eg, [log L(9\x, Y)] - Eg, [log k(Y\9, x)]. 

Let Q{9\9') = Eg,\logL(0\x,Y)}. It has been proved that 

log L(0\x) - ]agL(0®\x) > Q(9\0®) - Q(0«|0«) 

with equality only if 9 = 0®, or if k(Y\x, 0&) = k(Y\x, 9) for some other 9 ^ 9^. Choosing 
= argmaxgQ(9\9^) will make the difference positive and so, the likelihood could only 
increase in each step. The Expectation Maximization Algorithm is usually stated formally 
like this: 

E-step: Calculate function Q(9\0®). 

M-step: Maximize Q(0\0^) with respect to 9. 

The convergence of the EM sequence {9^} depends on the form of the likelihood L(9\x) 
and the expected likelihood Q{0\0^ 1 >) functions. The following condition, due to Wu ([21J), 
guarantees convergence to a stationary point. It is presented here as cited in [3]. 

If the expected complete-data log likelihood Eg/[logL(9\x, Y)\ is continuous both in 9 and 
9' , then all limit points of an EM sequence {9^} are stationary points of L(9\x), and L(9^\x) 
converges monotonically to L(0\x) for some stationary point 9. 

3.2 Derivation of an EM Algorithm for MTBP 

Let x be the observed set of particles, tt is the unobserved tree structure and p is the set of 
parameters - the offspring probabilities p(T v — > A) (the probability that a particle of type 
T v produces the set of particles ^4). Then the likelihood of the "complete" observation is: 

L(pK, x) = P(tt, x|p) = J] p(T v -> A) C(T ^ A ^ , 

v,A>Tv-*A 
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where c is a counting function - c(T v —> A; it, x) is the number of times a particle of type T v 
produces the set of particles A in the tree n, observing x. The probability of the "incomplete" 
observation is the marginal probability P(x\p) = P(w, x\p). For the EM algorithm we 
need to compute the function 

Q(p|p«) = E p(i) (log P(tt, x\p)) = P(*\x, P (0 ) Pi*, a?|p) 

7T 

= ^P(vr|x,p«) c(T„ ->• .4; 7T, a:) logp(T„ ->• .A) 

= £? p(i) c(T v -> A) log p(T v -> A). (2) 

X„-h4 

We need to maximize the function Q under the condition ^2^p(T v A) = 1, p(T v — > 
^4) > for every i; and A The Lagrange method requires to introduce the function 

$(p) = s P « c ( r « -> ^) logp{T v ->A) + A(l - £p(r v -> A)). 

T„-M. .4 

Taking partial derivatives with respect to p(T v —> A) and obtaining the Lagrangian multiplier 
A = ^2j[E d)(T v — > A), we get to the result that the re-estimating parameters are the 
normalized expected counts, which look like the ML estimators in the "complete" observation 
case Q, but the observed counts are replaced with their expectations. 

where the expected number of times a particle of type T v appears in the tree ir is: 

E p(i) c{T v ) = Y J P(K\x,p^)c{T v ;ir,x), 

and the expected number of times a particle of type T v gives offspring A in the tree ir is: 
E p(i) c{T v ^A) = Y P(*\x, P {i) HT v -> A; vr, x). 

IT 

It is easy to check that in this case the convergence condition stated above is fulfilled. We 
consider the representation 

Q(p|p (i) )= Y ^2 p ^\ x ^ p^) c ( t v A ^ x ) iogp(r v -> A), 

where P(ir\x, p^) is a polynomial function of pW-s - the offspring probabilities, estimated 
on the i-th step. Then Q(p|p^) is a sum of continuous functions in all the parameters p and 
pW, so it is also continuous. This way we are sure to converge to a stationary value, though 
it might not always be a global maximizer. 

It is a case of EM where the M-step is explicitly solved, so the computational effort will 
be on the E-step. The problem is that in general enumerating all possible trees tt is of 
exponential complexity. The method proposed below is aimed to reduce complexity. 
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3.3 The Recurrence Scheme 

We have shown that the E-step of the algorithm consists of determining the expected number 
of times a given type T v or a given production T v — > A appears in a tree it. A general method 
will be proposed here for computing these counts. The algorithm consists of three parts - 
calculation of the inner probabilities, the outer probabilities and EM re-estimation, which 
are shown below. 

Let us define the inner probability a(I,v) of a subtree rooted at particle of type T v to 
produce outcome I = {i\, 12, ■ ■ ■ , id} where it is the number of objects of type k (fig. 3). 
From the basic branching property of the process we get the following recurrence: 

a(I,v) = Y,P( T v^{T wl ,...,T Wk }) Yl a(l 1 ,T wl )...a(l k ,T Wk ) 

w Ii+...+Ifc=I 

where w = {T Wl , . . . ,T Wk } are all possible sets of particles that T v can produce. 



wl 'w2 • • • l w k 




Figure 3: The inner probabilities recurrence. 

The outer probability (3(1, v) is the probability of the entire tree excluded a subtree, 
rooted at particle of type T v and producing outcome I = 12, ■ ■ ■ , id} (fig- 4). The recur- 
rence here is: 

(3(i, v) = J2Y,p^ ^ T «( 2 ) > ■ ■ ■ > ^V) » 

W V 

X ^ @(I + J,W) ^ «(J2,W(2))"-a(Jm,V)) 
JCX— I J 2 +...+J m =J 

where {T v , v} = {T V ,T V ^, . . . ,T V{m) } are all possible sets of particles that T w can produce 
and X = {x\,X2, ■ ■ ■ , Xd} is the observed set of objects. 

For every T v G 7r the expected number of times c(T v ) that T v is used in the tree n could 
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Figure 4: The outer probabilities recurrence. 

be presented as follows: 

E p c(T v ) = p («\x, p)<T v ;tt, x) = ^ P ^'*^ c(T v ;ir, x) 

TV TV \^\P) 

= p&a lW ' t ' ) ' 

where a (I, t>) and are the inner and outer probabilities for all I C X. 

Similarly, the expected number of times a production T v — > {T Wl , . . . ,T Wk } is used could 
be calculated: 

E p (i)c{T v -)• {T Wl ,. . .,T Wk }) 

= J2 Yl P{l,v)a(I 1 ,w 1 )...a(I k ,w k )p^(T v -^{T wl ,...,T Wk }). 
I ii+...+i fc =i 

Dividing the expectations above, we obtain the EM re-estimation of the parameters: 
^ {T v {T Wl , . . . , T Wk Y) 

E E /8(I,«)a(Ii,«;i)...a(Ifc,«;fc)p«(r«->{T t01 ,...,r 1<;fc }) 
I Il+...+Ifc=I 

Ea(I,«)/3(I,«) 
i 

For several observed sets of objects the expected numbers in the nominator and denomi- 
nator are summed for all sets. 

Such generally stated the algorithm still has high complexity. In practical applications 
though, often there are small number of types and not all of the productions in the offspring 
distributions are allowed. Thus, for a number of specific dynamic programming 

schemes based on the above recurrence could be proposed, which will be less complex. 
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4 A Comprehensive Example 

We consider a MTBP with four types of particles - two nonterminal T±, T2 and two terminal 
T^, T 2 , and the productions that are allowed (with nonzero probability) are: 

tA^Tx}, T^{T u T 2 }, rA{7f}, T 2 %,T 2 }, T 2 i{Tj}. (4) 

To cover the possibility of observing particles of types T\ and T2 the additional productions 

Ti — K^i} and T2 — K^} are introduced. 

Now, let X = {1,0, 1, 1} and the process has started with one particle of type T\. We 
don't have further information about the distribution, so an uniform one is assumed for every 
type: p\ x = p\ 2 =p l T =p\ = 1/4 and p\ 2 =p\=p\ = 1/3. 

Let a({ii, 12, ij , i 2 }, v) be the probability of a tree rooted at a particle of type T v to 
produce the number of particles of types Ti,T 2 ,T± \T 2 respectively. The initialization of a 
should be as follows: 

«(iAO,o) =P\ = V 4 . "(0,0,1,0) =Pt = V 4 , "(0,1,0,0) = Pi = V3, "(0,0,0,1) = Pt = 1/3- 

For the level of two particles we are interested in all possible subsets of {1,0, 1, 1} con- 
taining two 1-s. The values of a are calculated below: 

"(1,0,1,0) = Pll"(l,0,0,0)"(0,0,l,0) + Pl2"(l,0,0,0) "(0,0,1,0) + Pl2 a (0,0,l,0)"(l,0,0,0) 

= 1/4.1/4.1/4 + 1/4.1/4.0 + 1/4.1/4.0 = 1/64, 

Similarly, 

"(1,0,0,1) = V48 and a( , ,i,i) = 1/48. 

Also 

"(1,0,1,0) = °> "(1,0,0,1) = °) and "(0,0,1,1) = °- 

Finally 

"(1,0,1,1) = Pll["(l,0,l,0)"(0,0,0,l) + "(1,0,0,1)"(0,0,1,0) + "(1,0,0,0)"(0,0,1,1)] + Pl2 ["(1,0,1,0) "(0,0,0,1) 

1 1 2 1 1 2 _i_ 2 1 _i_ 2 1 _L 2 1 1 

+ "(1,0,0,1)"(0,0,1,0) + "(1,0,0,0)"(0,0,1,1) + "(1,0,1,0)"(0,0,0,1) + "(1,0,0,1)"(0,0,1,0) + "(1,0,0,0)"(0,0,1,1)J 

= 1/4.(1/64.0 + 1/48.1/4 + 1/4.1/48) + 1/4.(1/64.1/3 + 1/48.0 + 1/4.0 + 0.0 + 0.1/4 + 0.1/48) 
= 1/256. 

Next follow the calculations of the outer probabilities (3({ii, 12, if, ijf}, v )- The initial 
values are 0, 1, 1}, 1) = 1 and /3({1, 0, 1, 1}, 2) = 0. Then: 

0(1,0,1,0) = /5(i,o,i,i)bn"(o,o,o,i) +Pi2«(o,o,o,i)] = 1-(V 4 .0 + 1/4.1/3) = 1/12, 
4,0,0,1) = 1/16, 4,0,1,1) = 1/16, 

4,o,i,o) = 4,o,i,i)^i2"(o,o,o,i) + 4,o,i,i)^2"(o,o,o,i)= 1-1/4.0 + 0.1/3.1/3 = 0, 
4,o,o,i) = 1/16, 4,o,i,d = 1/16, 
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0(1,0,0,0) - ^(1,0,1,0) [Pii a (o,o,i,o) + Pl2 a (0,0,l,0)] 

+ ^(1,0,0,1) [Pll Q (0,0,0,l) +Pl2«(0,0,0,l)l + 0(1,0,1,1) [Pil a (0,0,l,l) + ^20(0,0,1,1)] 

= 1/12.(1/4.1/4 + 1/4.0) + 1/16.(1/4.0 + 1/4.1/3) + 1.(1/4.1/48 + 1/4.0) = 1/64, 
0(0,0,1,0) = V64 and /^o,!) = 3/256, 

2 1 1 2 2 2 1 r 1 1 

^(1,0,0,0) = 0(1,0, l,0)Pl2 a (0,0,l,0) + 0(l,O,l,O)^22 a (O,O,l,O) + 0(1,0,0,1) bl2 a (0,0,0,l) 

+ 0fl,O,O,l)^22 a (O,O,O,l)] + 0(1,0,1,1) bl2 a (0,0,l,l) + 0(l,O,l,l)P22 a (O,O,l,l)] 

= 1/12.1/4.1/4 + 0.1/3.0 + 1/16.1/4.0 + 1/16.1/3.1/3 + 1.1/4.1/48 + 1.1/3.0 = 5/288, 
0(0,0,1,0) = 5/288, and /3 ( 2 OAOjl) = 1/256. 

Using that P(x\6) = q(X, 1) = a({l, 0, 1, 1}, 1) = 1/256, we are able to compute the 
expected values we need: 

EeciT,) —E^i 1 = ^(1/4-1/64 + 1/4.1/64 + 
+ 1/64.1/12 + 1/48.1/16 + 1/48.1/16 + 1/256) = = 4, 

E e c( Tl -+ { Tl ,T l} ) = p^)E i E ti«P( T i {ri.Tx}) = ^| = 1 
Eec( Tl -> { Tl ,T 2 }) = p^E i E ti^MATi -> {Ti,T 2 }) = ^| = 1 

^ ^ = TW)^ m -* {TT}) = vW 1/41/64 = if? = x > 
E ^ - ™ = pw/^ - ™ = 17k 1 / 41 / 64 = 172I = 1 

Thus, the estimated distribution for Ti is 

! _ £ eC (r 1 ->{r 1 ,r 1 » , _ e qC (t x -> {r^}) _ 

^mj _1/4 ' Pl2_ ^5 _1/4 ' 

which is the same as the initially chosen one, so this would be the final estimation. 
For the offspring distribution of T2 similar computations lead to the result: 

E e c(T 2 ) = l, E e c(T 2 -> {T 2 ,T 2 }) = 0, £ e c(T 2 -> {T 2 }) = 0, £7,c(T 2 -> {T 2 T }) = 1, 



so the estimation is: 



P22 = °> P2 = °> = !, 



which converges on the next iteration also. 



5 Simulation study 



Simulation experiment has been performed to study behaviour of the estimates obtained via 
the algorithm. Observations have been simulated according to the model Q in the previous 
section with offspring probabilities p\i = p\ 2 = v\ = 1/3 and p\ 2 = Pt = V^- Estimation 
has been performed using different sample sizes both for the number of observations, and the 
tree sizes as well. All the computations were made in R [19]. 

It is important to investigate how the size of the tree, which corresponds to the size of 
the "hidden" part of the observation, affects the estimates. In Table [l] are shown the result 
for small tree sizes and sample size 20. The most accurate estimates are obtained through 
averaging these results. It can be seen that there is great variation in the estimate of the 
individual distribution of type 2, thought the mean is close to the real values. For larger 
sample sizes the variance of the estimates is reduced, but there is some bias in the estimate 
for type 2 (Table [2]). Larger sample trees also lead to biased estimates for the individual 
distribution for type 2 (Table [3]). Using larger trees is also computationally more expensive. 



size 20 


Pt 


Pn 


P12 




Pt 


P 2 22 


size 20 


Pt 


Pn 


P12 


2 

Pt 


P22 


s.l 


0.35 


0.38 


0.27 


0.25 


0.75 


s.9 


0.34 


0.39 


0.27 


0.00 


1.00 


s.2 


0.26 


0.45 


0.30 


0.56 


0.44 


s.10 


0.36 


0.29 


0.36 


0.77 


0.23 


s.3 


0.33 


0.35 


0.32 


0.47 


0.53 


s.ll 


0.31 


0.29 


0.40 


0.46 


0.54 


s.4 


0.30 


0.33 


0.37 


0.47 


0.53 


s.12 


0.33 


0.41 


0.26 


0.59 


0.41 


s.5 


0.33 


0.38 


0.28 


0.00 


1.00 


s.13 


0.48 


0.17 


0.36 


0.92 


0.08 


s.6 


0.39 


0.42 


0.19 


0.00 


1.00 


s.14 


0.34 


0.25 


0.41 


0.64 


0.36 


s.7 


0.38 


0.35 


0.27 


0.94 


0.06 


s.15 


0.20 


0.33 


0.47 


0.93 


0.07 


s.8 


0.42 


0.25 


0.33 


0.72 


0.28 


s.16 


0.18 


0.25 


0.57 


0.66 


0.34 


mean 


0.34 


0.36 


0.29 


0.43 


0.57 


mean 


0.32 


0.30 


0.39 


0.62 


0.38 


st. dev. 


0.05 


0.06 


0.05 


0.33 


0.33 


st.dev. 


0.09 


0.08 


0.10 


0.30 


0.30 




all samples 


p\ 


P\\ 


P\2 


P T 


P22 




avg. 


0.33 


0.33 


0.34 


0.52 


0.48 


st.dev. 


0.07 


0.08 


0.09 


0.32 


0.32 



Table 1: Estimation obtained using small tree samples of size 20. 

The bias in the estimate for type 2 is due to the greater uncertainty in the process for type 
2: these particles could be generated by a particle of their own type, as well as, by a particle 
of type 1. For example, production of one particle of type 1 and two of type 2 could happen 
in two ways: once T\ — > {Ti,T2} and then T 2 — > {T2,T2}, or twice T\ — > \T\,T2\. So, in 
general, productions T\ — > {Ti,T2} take part more often in the estimation than productions 
T 2 — > {T2,T2}. As the branching is hidden and all possible generations have to be taken 
in account, this results in underestimation of p 22 an d some overestimation of p\ 2 when that 
hidden part gets larger. 
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size 50 


Pt 


Pn 


Pl2 


Pt 


P 2 22 














s.l 


0.29 


0.40 


0.31 


0.39 


0.61 














size 100 


Pt 


P\i 


Pl2 


Pt 


P22 


s.2 


0.34 


0.35 


0.31 


0.45 


0.55 


s.l 


0.31 


0.38 


0.31 


0.41 


0.59 


s.3 


0.36 


0.36 


0.27 


0.82 


0.18 


s.2 


0.38 


0.34 


0.28 


0.68 


0.32 


s.4 


0.39 


0.32 


0.30 


0.52 


0.48 


s.3 


0.29 


0.39 


0.31 


0.67 


0.33 


s.5 


0.31 


0.35 


0.34 


0.55 


0.45 


mean 


0.33 


0.37 


0.30 


0.59 


0.41 


s.6 


0.27 


0.43 


0.30 


0.81 


0.19 


st.dev. 


0.04 


0.03 


0.02 


0.16 


0.16 


mean 


0.33 


0.37 


0.30 


0.59 


0.41 














st. dev. 


0.04 


0.04 


0.02 


0.18 


0.18 














Table 2: Estimation obtained using small tree samples of size 50 and 100. 


size 50 


Pt 


Pu 


Pl2 


Pt 


P 2 22 














s.l 


0.33 


0.32 


0.35 


0.64 


0.36 


size 100 


Pt 




Pl2 


Pt 


P 2 22 


s.2 


0.33 


0.30 


0.36 


0.72 


0.28 


s.l 


0.34 


0.31 


0.36 


0.71 


0.29 


s.3 


0.33 


0.28 


0.39 


0.80 


0.20 


s.2 


0.33 


0.30 


0.36 


0.74 


0.26 


s.4 


0.34 


0.33 


0.33 


0.68 


0.32 


mean 


0.34 


0.31 


0.36 


0.73 


0.27 


mean 


0.33 


0.31 


0.36 


0.71 


0.29 


st.dev. 


0.00 


0.00 


0.00 


0.02 


0.02 


st.dev. 


0.00 


0.02 


0.02 


0.07 


0.07 





Table 3: Estimation obtained using larger tree samples of size 50 and 100. 



6 Conclusion 

A general EM algorithm has been proposed to find ML estimation of the offspring proba- 
bilities of MTBP with terminal types when only an observation of the generation at given 
moment is available. The example presented shows that the algorithm is straightforward and 
convenient to apply for a particular model. Simulation study shows that better estimates are 
obtained using smaller samples. Such algorithms would be useful in biological models like cell 
proliferation, genetics, genomics, evolution, and wherever a model of MTBP with terminal 
types is suitable. 
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